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Deriving ab initio model Hamiltonians for molecular crystals 


A. C. Jacko 

School of Mathematics and Physics, The University of Queensland, Brisbane, Queensland, 4072, Australia 

Developing realistic and precise models of the electronic properties of organic molecular crystals 
is crucial for understanding the full range of strongly correlated phases that they exhibit. By us¬ 
ing ab initio model construction methods, one can obtain unbiased non-interacting models of such 
systems from density functional theory, upon which one can base further (many-body) models. We 
will discuss the utility and advantages of ab initio model construction using Wannier orbitals. We 
will briefly review the approach, and then explain why it is so well suited to molecular crystals in 
particular. We discuss the ab initio construction of both non-interacting and interacting Hamitoni- 
ans, and highlight recent examples where such first principles models lead to importantly different 
results than fitted models. 
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BACKGROUND 
Molecular Crystals 

Molecular crystals, ordered periodic arrays of 
molecules, are known to exhibit a wide range of quan¬ 
tum mechanical phenomena, including unconventional 
superconductivity, quantum criticality, frustrated anti¬ 
ferromagnetism, and quantum spin liquid behaviour 
[8] . In some cases many of these phases can be found in a 
single material by tuning an external parameter, such as 
pressure or magnetic field. Often, one can control which 
phase is expressed by making subtle physical or chemical 


changes to the molecules mg. Along with the flexibility 
of the interactions within individual molecules, molecular 
crystals also have a range of intermolecular interactions. 
It is the subtle competition between the many intra- and 
inter-molecular interaction energies that brings the wide 
variety of phases seen in experiments so close together. 
These crystals tend to have a low effective dimension, and 
this likely contributes to the close competition between 
the phases imii]. 

Molecular crystals are an exciting testing ground for 
finding and understanding new emergent states of mat¬ 
ter. They often display competition between multiple 
emergent strongly correlated ground states OH). This, 
and the flexibility of organic chemistry, means that the 
emergent physics is often tuneable by subtle chemical and 
physical modifications (making slight variations around a 
core motif) [4]-[8|. In one notable case, a superconducting 
state in an organic molecular crystal is destroyed by sub¬ 
stituting some hydrogen atoms for deuterium, its heavy 
isotope m- That this extremely subtle change can have 
such profound consequences is both exciting and intim¬ 
idating. On one hand, it presents the inviting prospect 
of creating strongly correlated materials with technolog¬ 
ically desireable properties; on the other, the level of de¬ 
tail required to correctly predict the phase of a material 
can be substatial. 


Fig. [T] shows the phase diagram for a family of or¬ 
ganic crystals called Fabre salts [salts of TMTTF (shown 
in Fig. and an anion], which can be tuned through 
many different phases by applying physical pressure, or 
by slightly changing their anions, often thought of as a 
chemical pressure. As the size of the anion decreases, the 
TMTTF molecules pack closer together, as they would 
under the application of physical pressure. This range 
of accessible phases indicates that the many competing 
interactions are very close in energy. The many compet¬ 
ing energy scales in molecular crystals gives us access to 
phases with various interesting physical properties. 
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FIG. 1: Temperature-pressure phase diagram for the Fabre 
(TMTTF) aud Bechgaard (TMTSF) charge trausfer salts 
[3 E]? highlightiug the mauy accessible strougly-correlated 
phases. The phases showu are autiferromaguetic (AF), charge 
ordered (CO), Mott iusulatiug (MI), spiu Peierls (SP), spiu 
deusity wave (SDW), supercouductiug (SC) aud oue dimeu- 
sioual (ID), 2D, & 3D metals. The ambieut pressure positiou 
for each salt is iudicated with au arrow above the diagram. 



Predictive versus postdictive modeling 


phases seen in the Fabre salts. What one ideally would 
like is a systematic way of constructing an effective many- 
body Hamiltonian from first principles. Constructing the 
non-interacting part from first principles is currently em- 
minantly possible: By producing localised ‘Wannier’ or¬ 
bitals (discussed in more detail later), one can use the 
results of density functional theory (DFT) to construct 
a tight-binding lattice without first assuming its form. 
Here we give a brief overview of the theoretical develop¬ 
ment of the concept and practical details of using Wan- 
nier orbitals. For a much more complete and mathemat¬ 
ical account, turn to the excellent review of Marzari et 
al [12]. 

On a related note, it is worth commenting on the dis¬ 
tinction between ‘‘ab initio^ and ‘first principles’; ab ini¬ 
tio implies no empirical input, just calculations on the 
grounds of the many-electron Schrodinger equation us¬ 
ing the fundamental constants of nature such as Planck’s 
constant, the charge of the electron, etc. On the other 
hand, first principles allows for empirical parameters. 
Both density functional theory and Wannier orbital con¬ 
struction are in principle ab initio^ however particular 
implementations tend to include empirical parameterisa- 
tions (specific density functionals, for example) that are 
properly considered first principles rather than ab initio. 


Development of Wannier Orbitals 

In 1937 Gregory Wannier introduced the idea of con¬ 
structing localised sets of wavefunctions by fourier trans¬ 
forming Bloch states m- For a Bloch wavefunction for 
band n, '0n(k, r), the corresponding Wannier function for 
band n is 


There is an important philosophical point to be made 
here. The goal of science should be to make predictions 
about the nature of nature, and then to test those predic¬ 
tions. In the case of building models of materials, what is 
typically applied is a postdictive approach; one ‘knows’ 
that to have a model with the observed behaviour, it 
should be this lattice with that type of interaction (e.g. 
Heisenberg model on a triangular lattice, extended Hub¬ 
bard model on a square lattice, and so on). Thus, there is 
limited information gained about the system, whether 
one finds the behaviour one was searching for in such a 
model or not. One one hand, one chose the model phe¬ 
nomenologically, so finding the correct phenomenology is 
not profound. On the other, not finding the expected 
behaviour could be due to any number of reasons from 
the profound to the trite. 

Postdictive approaches can be useful, but one must 
go beyond them to gain a deeper understanding of the 
important commonalities and differences within a class 
of materials. For example, such an approach does not 
show much promise for describing all of the multitude of 


$n,R(r)=/ d=*ke-*>^-*^V’n(k,r), (1) 

JFBZ 

where R is any combination of the crystal lattice vectors 
with integer prefactors, ^R(r) is localised in the unit cell 
located at R, and the integral runs over the first Brillouin 
zone (FBZ). These new wavefunctions have the advan¬ 
tages of atomic orbitals (such as locality) while enforcing 
orthogonality. This allows one to treat localised exci¬ 
tations of individual electrons in metallic materials on 
the same footing as the ‘bulk’ electrons (the delocalised 
Bloch states). 

By the 1950’s these wavefunctions were widely known 
as Wannier functions, and of great use in understanding 
the physics of excitations in crystals. In 1953, George 
Koster introduced two new methods for defining Wannier 
functions without first having to solve the Schrodinger 
equation, and allowing one to use these orbitals to com¬ 
pute energy bands in crystals El- 

Walter Kohn put Wannier functions on a rigourous 
analytical grounding in 1959, showing that one can al¬ 
ways find a unique, real, symmetry-preserving, and ex- 
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ponentially localised functions for a given single band 
m- In this work he showed (although not in so many 
words) that Wannier orbitals were the ideal basis for 
the recently-developed tight-binding method (closely re¬ 
lated to the Hiickel method used in chemistry) [161 Hi] • 
(Kohn continued working on Wannier orbitals, and in the 
mid-90’s used the locality of Wannier functions, and the 
consequence that their interactions should decay expo¬ 
nentially, to propose a density functional theory method 
that scales linearly with the number of atoms naiis].) 
Jacques Des Cloizeaux further expanded the mathemat¬ 
ical grounding of Wannier functions, and identified what 
would later become known as the disentangling problem: 
if bands overlap, it is difficult to construct Wannier func¬ 
tions for just one of those bands (requiring one to ‘disen¬ 
tangle’ the target band from the other bands it crosses) 
m- This remains a general challenge in using Wannier 
functions to this day [12]. Due to the practical diffi¬ 
culty of the disentangling problem, and the extra inde¬ 
terminancy introduced in the disentangling proceedure. 
Wanner functions were not of signficant help in compu¬ 
tational electronic structure theory until the 90’s, when 
approaches based on density functional theory were in¬ 
troduced [21]. 

Wannier orbitals in Density Functional Theory 

The key breakthrough in the application of Wannier or¬ 
bitals occured when Nicola Marzari and David Vanderbilt 
formulated a generalised approach for generating maxi¬ 
mally localised Wannier functions for the case of multiple 
bands m- Not only that, they also described a numer¬ 
ical algorithm to produce such orbitals based on Bloch 
functions sampled on a mesh of points in /c-space, such 
as would be the output of a typical DFT code. This 
allowed for computations of Wannier orbitals in realis¬ 
tic, non-trivial cases. They also suggested the approach 
of using Wannier functions to construct effective model 
Hamiltonians for strongly correlated electron systems. 

A few years later, Marzari and Vanderbilt, along with 
Ivo Souza, extended their original approach to allow for 
entangled bands. Together they introduced an efficient 
disentangling methodology [22] that requires no addi¬ 
tional information over a usual Wannier construction, 
just one additional assumption. That assumption is that 
the ‘character’ of the Wannier orbitals (the contributions 
from particular basis functions) should vary as smoothly 
and slowly as possible; this is enforced via minimising the 
change in character across the Brillouin zone. 

These works laid the foundation for the wide-spread 
computation of Wannier orbitals in DFT codes. In 2008, 
the code wannierQO was released [23]. Developed by 
Arash Mostofi, Jonathan Yates, Young-Su Lee, along 
with Souza, Vanderbilt and Marzari, this code is now 
widely used, designed to interface with any DFT code to 


produce Wannier orbitals. It is now used in Wannier or¬ 
bital construction in FPLO [24], WIEN2k [25], Quantum 
ESPRESSO [26], ABINIT [27], and Eleur [28], to list just 
a few of the more popular DET codes for crystals. 


SEPARATION OF ENERGY SCALES IN 
MOLECULAR CRYSTALS 

Here we discuss the separation of energy scales that 
commonly occurs in molecular crystals and how this 
aids the construction of a minimal set of Wannier or¬ 
bitals. Molecular crystals tend to have a separation of 
energy scales in their non-interacting states, while there 
is competition amongst many possible strongly correlated 
ground states in the full many-body treatment. Despite 
the advances made in disentangling procedures, it re¬ 
mains a highly challenging task to produce high quality, 
reliable Wannier orbitals from entangled bands. This is 
particularly important if one is concerned about captur¬ 
ing the fine features of the electronic structure, which 
can have significant effects on the many-body state, as 
I will discuss explicitly in the case of crystals based on 
Pd(dmit) 2 . Molecular crystals provide an exciting play¬ 
ground for applying Wannier orbital based techniques to 
their greatest potential, because one can bypass the diffi¬ 
culty and ambiguity of disentangling procedures, one can 
determine the significance of the fine features of the elec¬ 
tronic structure in determining the rich phase diagrams 
of these materials. 

The separation of energy scales in molecular crys¬ 
tals is straight-forward to understand: the molecules 
are held together by (strong) covalent bonds, while the 
crystal is held together by much weaker intermolecular 
forces; van der Waals, tt- stacking, and hydrogen bond¬ 
ing. The strong forces within a molecule produce a set of 
well spaced molecular orbitals (MOs), and these orbitals 
weakly couple between molecules, as illustrated in Fig. 

producing bands that are narrow on the scale of the 
MO energy gaps. 

We can understand this more concretely by consider¬ 
ing a toy example: a ID, two orbital tight-binding model. 
For simplicity we consider a chain of spinless fermions, 
with orbitals a and b on each site i, with nearest neigh¬ 
bour interactions. The Hamiltonian is 

H = ^^{na,i - nb,i) + Y (2) 

i a=a,b 

where A is the energy difference between the orbitals, 
and ta is the inter-site hopping for orbital a. In the case 
of molecular crystals, the orbitals are molecular orbitals 
of single molecules. The energy difference between the 
molecular orbitals (A) comes from the inter-atomic hop¬ 
ping within a molecule, typically a 7r-type overlap. A 
typical energy scale for this difference between molecular 
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FIG. 2: Discrete energy levels of a single molecule (left) (here, 
TMTTF sans hydrogens), and continuum band states of the 
resulting molecular crystal (right). The interatomic interac¬ 
tions within the molecule set the energy scale for the molec¬ 
ular orbitals, while the inter-molecular coupling determines 
the width of the bands. In typical cases these energy scale 
are quite different (as illustrated). 


orbitals in an organic molecule is a few eV (see for exam- 
pie [29]). The inter-site hopping comes from the overlap 
of molecular orbitals on different molecules, and is ex¬ 
ponentially suppresed by distance. These energy scales 
are on the order of 10 - 100 meV for nearest neighbour 
overlaps (see for example H)- Thus it is often the case 
in such systems that |A| > |ta| + |4|; the bands result¬ 
ing from each molecular orbital are narrow enough and 
well-separated enough that they do not overlap in energy. 
Thus, depending on filling, one can consider just one or¬ 
bital or the other as the foundation for an effective model 
Hamiltonian. 

In applying the Wannier construction procedure out¬ 
lined above, we have glossed over the details of limit¬ 
ing the Fourier transform window to some small energy 
range. In the simplest case of a single band system, it 
is clear that this originally-infinite window can be trun¬ 
cated to be exactly the bandwidth of the single band 
without any loss of generality. However, in multi-band 
systems, extracting a subset of bands can become diffi¬ 
cult. To understand why, let us consider again the ID, 
two orbital model. 

For |A| < \ta \ + 1^51, there are gapless excitations pos¬ 
sible between these bands; the two bands have weight in 
an overlapping energy range. Thus, even in this simplest 
case without band crossings or interactions, the presence 
of a second band makes picking out the states of the first 
band non-trivial. When these bands cross or hybridise. 



FIG. 3: Electronic properties of (TMTTF) 2 PF 6 at T = 4 K, 
reporduced from [7]. a) Band structure and density of states, 
b) path through /c-space, c) Fermi surface in the kz = 0 plane. 
The total density of states is shown with the solid black line, 
and the partial density of states of the anions (x 100) is shown 
in dashed orange. The partial density of states shows that the 
two bands at the Fermi level are nearly purely TMTTF, with 
large energy gaps on either side, demonstrating the separation 
of energy scales. 


this further complicates the procedure. This problem is 
very difficult to solve in general and is known as the dis¬ 
entangling problem. 

Now, the separation of energy scales in molecular crys¬ 
tals comes in to play. As discussed above, because of 
the often quite different energy scales of inter- and intra¬ 
molecular interactions, it is typical to find well isolated 
sets of bands in the band structure of a molecular crys¬ 
tal, as illustrated in Fig. This property means that 
one can bypass the difficulty and ambiguity of projective 
disentangling procedures. Thus minimal input is needed 
into the WO construction procedure in molecular crys¬ 
tals; one just inputs how many orbitals you would like, 
spanning what energy window. 


AB INITIO MODEL CONSTRUCTION 


It is important when modeling these systems that the 
models we use be as accurate and unbiased as possible; 
starting with preconceptions of how the model ‘should 
look can limit what one finds. 































































5 


Tight-Binding Models 

DFT can give us useful information about the non¬ 
interacting electronic properties of a system. To utilise 
this information, we will construct a tight-binding model 
from the DFT using a rigorous Wannier orbital construc¬ 
tion technique. This procedure creates localised orbitals 
that accurately represent the electronic properties found 
by DFT for the frontier electrons, those that determine 
the low-temperature physics. As discussed above, the 
separation of energy scales makes Wannier orbital con¬ 
struction straightforward in molecular crystals. Once one 
has local Wannier orbitals one can construct a first prin¬ 
ciples tight-binding model. 

First principles versus Fitting 

One might ask why is all of this effort justified? Why 
not simply write down a perfectly good tight binding 
model and fit it to a first principles band structure? (As 
is often done, see [8l l3Qti32] for just a sample.) For very 
simple systems, where the relevant tight-binding model 
is clear, and only has a few parameters, a first principles 
method is probably not justified as the band structure is 
well reproduced by fitting methods [8l [3QH32] . However, 
in the much more common situation of a somewhat am¬ 
biguous tight-binding model with an unknown number of 
relevant parameters, a first principles approach is ideal, 
as I will discuss further below. In fact, in some systems 
where a quite simple model seemed obvious, it has been 
shown that a somewhat more nuanced model does a bet¬ 
ter job of capturing the important many-body physics of 
the material (as discussed in the context of Pd(dmit) 2 , 
below). 

John von Neumann is famously claimed to have said 
“With four parameters I can fit an elephant, and with 
five I can make him wiggle his trunk.” [33] (this is more- 
or-less true, as shown in Fig. after Ref. [34|); this cap¬ 
tures the essence of the problem of fitting energy bands 
to a tight-binding model. With enough parameters in the 
fit, it is difficult to have a bad fit of the band structure; 
a model with many parameters might reproduce the dis¬ 
persion without having any connection to a realistic mi¬ 
croscopic description of the system. As such, a good fit 
provides almost no information, especially not about the 
quality of the microscopic model to which you are fitting. 
Worse, when there are many parameters, very different 
values can produce similarly good fits by whatever op¬ 
timisation metric you are using. Often, these different 
sets of values have importantly different physical con¬ 
sequences (for example changing the electronic dimeri- 
sation of a chain, or the localisation of charge). These 
differences can lead to importantly different many-body 
ground states, as will be discussed further below. 

To demonstrate this, I will discuss a particular case 



FIG. 4: A four-parameter fit to an elephant, produced as 
described in [34]. A fifth parameter does indeed allow it to 
wiggle its trunk. While conforming to the letter of the state¬ 
ment, this implementation somewhat defies its spirit as the 
parameters are complex numbers, thereby carrying twice the 
information of four real numbers. 


of producing a tight-binding model for an organic 
molecular crystal by fitting and via Wannier orbitals. 
TMTTF 2 ASF 6 is one of the Fabre salts, a family of or¬ 
ganic charge-transfer salts with a rich phase diagram, 
given in Fig. These fiat organic molecules tt- stack into 
one dimensional chains along the crystallographic a di¬ 
rection (shown in Fig. |^, with some inter-chain coupling 
in the a — b plane, and minimal coupling in the c direction 
(where the anionic AsFe layers introduce a large spacing 
between TMTTF layers). Thus this system is largely ID 
electronically, with some 2D nature introduced by next- 
nearest neighbour hopping and further terms 0171135]. 

In this family of materials, once one decides to in¬ 
clude any 2D hopping terms, one encounters the problem 
that there are many terms of similar magnitude; to avoid 
neglecting terms of the order one keeps, one needs to 
add many parameters to the tight-binding model. These 
many degrees of freedom cause problems for fitting pro¬ 
cedures; one finds many local minima with similar op- 
timiziation functions values, but very different parame¬ 
ters, with important physical consequences. This is the 
heart of the problem; that such fits are examples of sloppy 
models [36] : changes in one parameter can be almost com¬ 
pletely masked by compensatory changes in other param¬ 
eters. 

Figure 0 shows the results of 5000 runs of a fitting pro¬ 
cedure applied to the band structure of TMTTF 2 ASF 6 
[a El]; a pre-defined tight binding lattice is input, and 
the values of the 8 different hopping integrals (c.f. Fig. 
1^ are optimised to provide the best fit to the band struc- 
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FIG. 5: Crystal structure of TMTTF 2 ASF 6 at room temper¬ 
ature, showing a 7r-stacked chain along the a direction, and 
spacing due to anions in the c direction. Reproduced from 

m- 


ture (quantified by the least-squares error over the set of 
points the bands are sampled on). Each run starts with 
a (different) randomised set of fitting parameters and it¬ 
erated on. Due to the nature of this method, each of 
the 5000 fits is basically unique. Rounding each t to the 
nearest 0.1 meV, there are 26 unique fits; of the 5000 re¬ 
sults, 74 have the minimum objective function. A given 
run has a 1.5% chance (1/67) of finding this ‘best’ so¬ 
lution! There are many more fits with a just slightly 
higher value of the objective function, and these fits con¬ 
tain contradictory physical information. For example, 
these different fits make different predictions about the 
electronic dimerisation of the system; whether the elec¬ 
tronic dimers are on the structural dimers or not. This is 
an important difference, and can change as a function of 
pressure for example [7]. While most of the parameters 
have positive and negative equivalents, a few parameters 
are very precisely determined, and with a fixed phase. 
This means that the relative phases of the hopping in¬ 
tegrals cannot be absorbed by a gauge transformation; 
these different sets of parameters have different physical 
meanings. The number of parameters used in the fit also 
effects the outcome. Removing some t parameters from 
the inputed model will cause the other t values to change. 
Damningly, in this system the optimal tight-binding fit is 
quite different to the set of hopping integrals found from 
Wannier orbitals, as shown in Fig. The Wannier and 


fitted parameters make contraditary predictions about 
the electronic dimerisation of the system (as seen in the 
magnitude of the first two t’s). In the limit of including 
all Wannier overlaps out to infinite distance, the set of 
Wannier tight-binding parameters will exactly reproduce 
the band structure when the bands are not entagled. 

Overall, it is hard to take fitting procedures too seri¬ 
ously with more than a couple of parameters in the fit. 
With small numbers of parameters, fitting becomes more 
stable, and in certain systems a few parameters is enough 
to acurately describe the band structure [H [30l |32] . Even 
then, a fitted model should be considered an effective 
model that has potentially lost important detail in ‘in¬ 
tegrating out’ the full set of parameters. In a sense 
these are ‘variational Hamiltonians’; they are optimised 
by some metric, but there is no assurance that they rep¬ 
resent the underlying microscopic physics. 

By producing Wannier orbitals for molecular crystals, 
and computing a set of t’s from those, one finds a sin¬ 
gle set of parameters that is reliable and robust; one 
can believe them just as much as one believes the other 
results of the DET computation. In the general case 
of Wannier orbital construction, there is ambiguity in¬ 
volved in disentangling bands to produce the desired set 
of Wanniers m- When the bands one cares about are 
well-separated from the bulk, this ambiguity is gone. 
Not only that, but one can gain knowledge by look¬ 
ing at the Wannier orbitals themselves. In the case of 
TMTTE 2 ASE 6 , the Wannier orbitals are localised to sin¬ 
gle TMTTE molecules, as shown in Eig. Here, the 
Wannier orbital is very much like the HOMO (highest 
occupied molecular orbital) of a single TMTTE molecule 
in vacuum. Having this real-space orbital allows us to 
do many further computations based on the DET, as 
well as producing a single robust parameter set for a 
tight-binding model (Eig. [^. This Wannier based model 
construction technique is being used more and more in 
molecular systems [illllslHio]. 


Including Interactions 

Correctly parameterising many-body effects is an im¬ 
portant and challenging task; it is the competition be¬ 
tween energy scales that leads to the interesting physics 
in many systems, so small relative changes in large pa¬ 
rameters can have large effects mi- Often, these param¬ 
eters are estimated without careful consideration of the 
assumptions involved. Eor instance, if one considers a 
Hubbard model on a dimer with an inter-monomer hop¬ 
ping t, on-monomer Hubbard repulsion and inter¬ 
monomer Hubbard repulsion V; in the limit Um 00 , 
E ^ 0, the effective Hubbard repulsion in the dimer 
orbitals is Ud = 2t [42j [43]. This assumption is often 
used, since it allows one to estimate many-body param¬ 
eters from straightforward band structure calculations. 
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FIG. 6: Instability of fitting with many parameters. 5000 
fits of an 8 t tight-binding model to the band structure of 
TMTTF 2 ASF 6 produced in [7]. Each line is a set of parame¬ 
ters resulting from one run of the fitting algorithm. The best 
fit parameters are shown in red. They are inconsistent with 
the parameters found from Wannier orbital overlaps (given in 
green). Only 1.5% of the runs found the minimal value of the 
objective function. The histogram shows the sets of minimi¬ 
sation function values produced in the 5000 runs. Less than 
6% of runs are in the ‘best’ segment of the histogram, and 
1/4 in the best two segments. 


or molecular Hiickel calculations. However, it is not well 
justified. It has since been shown that in the more realis¬ 
tic case of Um ^ H ^ t, that Ud = [44]. Thus, 

one still needs to be able to correctly compute many-body 
parameters to estimate the dimer parameters (even be¬ 
fore considering screening). None-the-less, this approxi¬ 
mation continues to be used (for example [32l 05] 1^). 
often without stating the strong underlying assumptions. 



FIG. 7: Wannier orbital for TMTTF 2 ASF 6 . Note that this 
orbital is localised to a single molecules, and very much like 
the HOMO of an isolated TMTTF molecule. Having this 
real-space orbital allows us to do many further computations 
based on the DFT, such as computing a tight-binding model 
by taking real space overlaps of such orbitals. Reproduced 
from [Tj. 


One might think that, given these real-space Wannier 
orbitals for a particular system, it must be straightfor¬ 
ward to calculate the many-body Coulomb integrals and 
parameterise a Hubbard model. However, computing 
these terms by simply evaluating the Coloumb energy 
for each orbital neglects screening (equivalently, relax¬ 
ation of the bulk states). Screening can easily suppress 
the Hubbard U by an order of magnitude [47] . 

In classical electromagnetics, there are many tech¬ 
niques for determining the response of a bulk to a per¬ 
turbing field (analogous to the case here of computing 
the screening/relaxation of a doubly occupied orbital). 
The discrete dipole approximations (also called the cou¬ 
pled dipole approximation) is one such technique. In this 
approximation, one discretises the bulk as a set of polaris- 
able dipoles, and self-consistently solves their response to 
the perturbing field and to each other [48| . This method 
is very much like a technique applied to molecular crys¬ 
tals to compute screened Coloumb parameters. By rep¬ 
resenting each molecule by a set of polarisable (classical) 
dipoles, and placing a perturbing charge on one lattice 
site, one can compute the correction to the Hubbard re- 
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FIG. 8: Tight-binding lattice for TMTTF 2 ASF 6 produced 
from Wannier orbitals. The thickness of the lines is propor¬ 
tional to the magnitude of the t; to = 175 meV, ti = 157 
meV, and the rest are < 25 meV. Reproduced from [7]. 


pulsion due to the polarisation of the rest of the molecules 
in the crystal [49]. This technique, though promising, 
has only been applied to a single molecular crystal (TTF- 
TCNQ), with no new applications apparent in the 5 years 
since the original publication. 

An alternative approach to computing screened 
Coulomb parameters from first principles has gained 
prominence recently, the constrained random phase ap¬ 
proximation (cRPA) [50]. This technique is also based 
around computing the polarisation of the system, but in 
this case, the quantum mechanical polarisation function 
in the random phase approximation. Here we discuss 
RPA, cRPA, and its application to molecular crystals in 
more detail. Practically, it also relies on having Wannier 
orbitals for a few relevant bands, and so like the tight- 
binding model construction it is particularly suitable to 
molecular crystals. 

Constrained Random Phase Approximation 

The random phase approximation (RPA) was intro¬ 
duced by David Bohm and David Pines in the 1950’s to 
include the effects of screening into models of electron 
gases [SHai. Murray Gell-Mann and Keith Brueckner 
placed this approximation on a firmer footing, showing 
that the RPA can be derived from a self-consistent series 



+ 


FIG. 9: Random phase approximation bubble diagrams ap¬ 
propriate for calculating the polariation function. Bold lines 
are fully interacting Greens functions, while non-bold are non¬ 
interacting [56] . 


of leading order Feynman diagrams m, an example of 
which is illustrated in Fig. [^ 

Given a basis of occupied and unoccupied states, one 
can compute an RPA polarisation function 

occ unocc 

* J (3) 

( _ ^ ^ _ ^ _ _) , 

\UJ — 6j £i UJ €j — Si — J 

where i (j) runs over the occupied (unoccupied) single 
particle states. With this polarisation function one can 
compute the screening effects of this system. This is ex¬ 
actly what one wants if computing the effects of an im¬ 
purity in such a system, for example. However, if one 
wants to include many-body effects in a lattice model, 
then this constitutes an overcounting of the effect. 

In 2004, Aryasetiawan et al. introduced a new, pre¬ 
cise method for constructing sets of screened effective 
model parameters for strongly correlated lattice models 
m- The constrained random phase approximation (con¬ 
strained RPA or cRPA) is a systematic way of account¬ 
ing for screening in the many-body parameters computed 
for some basis orbitals. The system is divided into two 
fragments; the active subspace (often labelled ‘d’), the 
space spanned by the orbitals of interest; and the rest 
of the bands (labelled ‘r’). On a conceptual level, this 
procedure computes the effects of transitions involving 
the T’ subspace with RPA (by computing the polarisa¬ 
tion function due to these transitions), while leaving the 
transitions within the ‘d’ subspace to be dealt with in 
the many-body model that results m- This proceedure 
allows one to generate all the terms resulting from the 
Coulomb interaction, on- and off-site repulsive and mag¬ 
netic interactions. Practically, one constrains the sums in 
Eq. [^to exclude transitions within the active subspace. 

The partitioning idea at the core of cRPA works best 
in the same situation that the Wannier orbital procee- 
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dure itself works best: a set of relevant bands well sepa¬ 
rated from the bulk. In the situation of entangled bands; 
where the natural basis one would like to use mixes with 
bands due to other states; one can apply the Wannier dis¬ 
entangling proceedure to construct a disentangled basis 
m- In inorganic system, there are difficulties in disen¬ 
tangling the target bands from the bulk. Nonetheless this 
approach was quickly applied to transition metal systems 
and simple transition metal oxides miEzHsn]. 

This approach has been applied to only a few organic 
crystals [38l |39]. In those cases where it has, it finds 
sometimes importantly different parameter values. In the 
ET charge transfer salt /t:-(ET )2 Cu 2 (CN) 3 , cRPA pre¬ 
dicts a value of Ujt for the dimer about twice as large as 
was estimated from a Hiickel analysis of a dimer (using 
an optial conductivity estimate of the monomer value, 
Um), ^ 15 vs ^ y | 3 gj gj] ^ ^ 

simple Hubbard model, this would place this material 
well into the insulating phase, contrary to the observed 
metallic behaviour. The cRPA analysis also showed that 
off-site V terms are significant, V/U ^0.5, meaning that 
to properly understand the system one must consider an 
extended Hubbard model |38] . While the optical conduc¬ 
tivity estimate for the monomer Um is quite reliable, the 
assumptions in using this to estimate Ud are not. This 
Wannier-based approach gives us a reliable first princi¬ 
ples estimate of all the Hamiltonian parameters on the 
same footing. 

FIRST PRICIPLES APPROACH FINDS 
IMPORTANT DIFFERENCES 

Here we discuss particular examples to show that using 
a first principles approach can give importantly different 
results and insights than a standard fitting approach; be 
it caused by subtleties of parameter variations or quali¬ 
tatively different lattices. 

EtMe3Sb[Pd(dmit)2]2: Fine details matter 

To demonstrate the importance of finding a robust 
set of model parameters we will turn to the example of 
EtMe 3 Sb[Pd(dmit) 2]2 5 a spin-liquid candidate material 
and part of a family of organic molecular crystals with 
a rich phase diagram; as well as the spin-liquid phase, 
these materials have Mott insulating, superconducting, 
spin density wave and valence bond solid phases [2l|60^ 
[63] . Constructing a coherent picture of this family of ma¬ 
terials and their many phases is highly challenging. This 
effort has been hindered by the fact that, in the usual de¬ 
velopment of microscopic models, many approximations 
are made without fully understanding their consequences 
0 . 

The typical approach in EtMe 3 Sb[Pd(dmit) 2]2 and the 


related family of materials is to focus on a dimer model, 
where the dimers of Pd(dmit )2 sit on a t — t' triangu¬ 
lar lattice. Parameters are either fit or mapped to this 
non-interacting model before many-body effects are con¬ 
sidered (see for example Refs. [3TJ [63]). It has since 
been shown that a fully-anisotropic triangular lattice 
(FATE; t — t' — t") better represents the electronic struc¬ 
ture m HU El]. Further, it was shown that a FATE 
allows one to reproduce the observed many-body prop¬ 
erties, predicting a spin-liquid ground state for reason¬ 
able parameter values in EtMe 3 Sb[Pd(dmit) 2 ] 2 , while 
the t — t' model does not [3]. Fig. shows the 
phase diagram for the Hubbard model (as a function 
of U/t) on the isotropic triangular lattice, t — t' tri¬ 
angular lattice, and fully anisotropic triangular lattice 
(FATE), each with tight-binding parameters consistent 
with EtMe 3 Sb[Pd(dmit) 2]2 (computed with variational 
quantum Monte Carlo). First principles estimates pre¬ 
dict IJ It max ^11 [39]. The FATE enters the spin-liquid 
phase at this point, while the t — t' and isotropic lat¬ 
tices would predict an insulating phase, with this value 
of U very far from the critical value. Generally, the extra 
anisotropy seems to destablise the insulating phase rel¬ 
ative to the metallic and spin-liquid phases. It is worth 
noting that these variational quantum Monte Carlo re¬ 
sults are not definitive; however, if nothing else, they are 
indicitive of the important consiquences that even slight 
parameter changes can have. 

Such highly anisotropic models have since become in¬ 
creasingly used in investigations of this family of mate¬ 
rials [8]. Having reliable and believable predictions of 
the degrees of anisotropy in these materials (were the 
fine variation in parameter values can be attributed to 
physics and not a quirk of the particular fit one is ap¬ 
plying) will be vital for building an understanding of the 
whole class of materials. 


/^-(BEDT-TTF)2 salts: Long range terms 

The BEDT-TTF (bis(ethylenedithio)- 

tetrathiafulvalene, or ET) organic charge transfer 
salts are a family of quasi-2D crystals that exhibit 
a wide range of strongly correlated phases (such as 
non-BCS d-wave superconductivity) [U El (651168] . Un¬ 
derstanding this wide range of phases requires a good 
effective model and good model parameters. It was in 
these materials that the shortcomings of the U ^ Untra 
approximation (discussed above) were made clear, 
showing that it leads to a systematic underestimate of 
U [44]. Once a realistically large value of U is used 
(computed with cRPA and found to be a 50% - 100% 
increase over previous estimates), a straightforward 
Hubbard model of the dimer lattice does not capture 
anything but the Mott insulating phase [38|. These 
cRPA parameter estimates also showed that the nearest 
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Anisotropy 



FIG. 10: Phase diagrams of the Hubbard model on the 
isotropic triangular lattice, t — t' triangular lattice, and fully 
anisotropic triangular lattice (FATL), for parameters consis¬ 
tent with EtMe 3 Sb[Pd(dmit) 2 ] 2 . Note that the transition to 
the spin-liquid phase occurs for a much smaller value of U in 
the FATL. Phase diagram determined with variational quan¬ 
tum Monte Carlo [3]. Reproduced from [3]. 



FIG. 11: Lattices for Mo 3 S 7 (dmit) 3 : the phenomenological 
‘triangular necklace’ lattice on the left, and hrst priciples 
kagomene on the right. While quite different, both have in¬ 
teresting topological properties and can provide insights into 
the behaviour of Mo 3 S 7 (dmit) 3 . 


neighbour inter-site Coulomb interactions are significant 
{V/U ^ 0.5), and that they decay slowly with distance 

m- 

By applying first principles model building techniques, 
it was found that describing the phases of the ET salts 
requires models like the extended Hubbard model with 
significant and long-ranged inter-site interactions. This 
kind of model, although it has more parameters, it has 
no more free parameters. Additionally, the inclusion of 
long-range Coulomb terms has important implications for 
the energetics of ordered phases [38] . 

Mo3S7(dmit)3: An unexpected lattice 

In the previous examples, we showed how details in 
the model parameters are found to have significant con¬ 
sequences on the predicted ground state. We now turn to 
a system where, by using a first principles method, one 
finds a totally different lattice model than any previously 
considered for this system. 

Mo 3 S 7 (dmit )3 is a single component molecular crys¬ 
tal that was designed to be metallic. However, it was 
found to be an activated insulator with an activation en¬ 
ergy of 34 meV [69]. Further, it was found to have no 
sign of any magnetic order down to very low tempera¬ 
tures {J/ksT 50 [69|); a possible sign of a spin-liquid 
state. Based on the apparent ID physicical properties, its 
crystal structure, and initial bandstructure calculations, 
Mo 3 S 7 (dmit )3 was modeled with a one-dimensional lat¬ 
tice [69U7T] . This is the ID ‘triangular necklace’ lattice, 
illustrated in Fig. El The ground state of the Hub¬ 
bard model on this lattice at 2/3 filling (as appropriate 
for Mo 3 S 7 (dmit) 3 ) is found to be in the Haldane phase. 


consistant with the experimental evidence [Toiin]. 

However, Wannier orbital tight-binding model con¬ 
struction based on density functional theory for 
Mo 3 S 7 (dmit )3 predicts that at the single electron level, 
this system is actually 2D with coupling between the 2D 
layers. The lattice of the 2D layers is an unusual dec¬ 
orated honeycomb lattice, the ‘kagomene’ lattice (illus¬ 
trated in Fig. [IT] ); interpolating between the graphene 
(honeycomb) and kagome lattices [72]. These lattices 
have quite different properties, and provide quite dif¬ 
ferent pictures of the physics of Mo 3 S 7 (dmit) 3 . The 
kagomene lattice has been studied theoretically before 
[73H77], but never seen in a real system. This first 
priciples approach found a layered kagomene lattice in 
Mo 3 S 7 (dmit )3 quite unexpectedly, demonstrating the 
novel insights this appoach can yield. The microscopic 
picture produced is quite different from the phenomeno¬ 
logical model. 

The one dimensional behaviour can be understood on 
the grounds of the kagomene lattice: just like kagome, 
this lattice has exactly localised states m, illustrated 
in Fig. [^ Once the 2D kagomene lattice is extended 
into 3D, these localised states become ID bands. These 
emergent ID states are topological; their degeneracy de¬ 
pends on the boundary conditions of the lattice. Thus, 
dispite the hopping integrals having similar magnitudes 
in every direction (in fact, slightly smaller in the stack¬ 
ing direction), one recovers the ID behaviour and gains 
some important insights about the potential topological 
properties of this system. 

As a phenomenological model, the necklace lattice does 
a good job of reproducing the observed magnetic prop- 
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FIG. 12: Localised states on the kagomene lattice, showing 
two plaquette states, the anti-bonding Ar and bonding Br/ , 
and a topologically non-trivial loop state, which can only exist 
with periodic boundary conditions. Reproduced from m- 


erties of Mo 3 S 7 (dmit )3 [7T]. On the other hand, the 
kagomene model provides a natural explanation for the 
quasi-ID behaviour, and highlights the interesting topo¬ 
logical flat bands analogous to those seen in the kagome 
lattice m- In addition, one can find a lattice closely 
related to the necklace model as a limiting behaviour of 
an interacting model in the kagomene lattice, and the 
many-body behaviour of this model is very similar to the 
necklace model m- One can naturally find new terms 
to extend the necklace model in a consistent way by in¬ 
troducing higher-order terms in these limits, for example 
including the chiral next-nearest neighbour terms m- 

SUMMARY 

Over the last decade Wannier orbitals have become an 
important tool for predictive physics. By constructing 
Wannier orbitals for frontier bands, we can derive effec¬ 
tive models that avoid our biases. These models are ro¬ 
bust and reliable, and allow us to make detailed compar¬ 
isons between materials, and start to extract some gen¬ 
eral behaviours. They can also find models that we might 
never have expected to see, leading us to new insights. 
By moving to this kind of assumption-free methodology 
for model construction, we can move to a truly predictive 
approach. We can avoid the dangers of relying on vari¬ 
ational Hamiltonians, and allow ourselves to find truely 
unexpected things. 
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